Quantum Monte Carlo Simulation of the High-Pressure Molecular- Atomic Crossover 

in Fluid Hydrogen 
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A first-order liquid-liquid phase transition in high-pressure hydrogen between molecular and 
atomic fluid phases has been predicted in computer simulations using ab initio molecular dynam- 
ics approaches. However, experiments indicate that molecular dissociation may occur through a 
continuous crossover rather than a first-order transition. Here we study the nature of molecular 
dissociation in fluid hydrogen using an alternative simulation technique in which electronic corre- 
lation is computed within quantum Monte Carlo, the so-called Coupled Electron Ion Monte Carlo 
(CEIMC) method. We find no evidence for a first-order liquid-liquid phase transition. 

PACS numbers: 61.20.Ja, 62.50.+p, 64.70.Ja, 77.84.Bw 



The behavior of hydrogen under a range of thermo- 
dynamic conditions is an important problem in theoret- 
ical physics, as well as in planetary science and high- 
pressure physics. Hydrogen exhibits a rich variety of 
properties including several molecular and non-molecular 
phases, fH a fluid metal- insulator transition, and pos- 
sible fluid 3 and superconducting [J| phases at low tem- 
perature. There exists much uncertainty about the equi- 
librium properties of high-pressure phases, including the 
structure and relative stabilities of a number of solid 
phases, [1, @] the specific shape of the molecular-solid 
melting curve and the reason for the density maximum, [sj 
and the conditions of molecular dissociation in the fluid. 

Shock-wave experiments Q have hinted that the met- 
alization of hydrogen in the liquid state occurs in a con- 
tinuous manner at temperatures above 3000K. However, 
recent computer simulations [sl, 01 using Car-Parrinello 
molecular dynamics (CPMD) with density functional 
theory within the local density approximation (DFT- 
LDA) and generalized gradient approximation (DFT- 
GGA) predict that the molecular dissociation in the 
fluid happens through a first-order phase transition, with 
discontinuous metalization, at temperatures lower than 
3000K. These results could be reconciled by having the 
critical point for the transition at T < 3000K. In addi- 
tion, other simulations of hydrogen fluid phases, which 
do not directly address the nature of the molecular dis- 
sociation, have been reported 

In this Letter, we present the results of a study 
of this transition employing a new quantum Monte 
Carlo (QMC) method, coupled electron-ion Monte Carlo 
(CEIMC). The first QMC approaches for studying 
high pressure hydrogen were limited to studying solid 
phases at T = OK: variational Monte Carlo (VMC) 
and diffusion quantum Monte Carlo (DMC)[ll|, IT^ . 
Finite-temperature restricted path integral Monte Carlo 
(RPIMC) was used to investigate the possibility of a hy- 
drogen liquid-liquid phase transition [13|, finding a con- 
tinuous molecular dissociation at T ~ lOOOOK in agree- 



ment with experiments. However, RPIMC becomes in- 
efficient and inaccurate at lower temperatures. 

The CEIMC approach, in contrast with earlier QMC 
methods, involves simultaneous evolution of separate but 
coupled Monte Carlo simulations for the electronic and 
ionic subsystems, which may then be treated on differ- 
ent levels of approximation. A large gain in computa- 
tional efficiency can be made for lower temperatures by 
employing the Born-Oppenheimer approximation (BOA) 
and requiring that the electronic subsystem remain in its 
ground state (T — OK). The ionic system is treated at 
finite temperature either as a set of classical particles, 
or quantum mechanically using the imaginary-time path 
integral formalism. This decoupling is good provided the 
thermal excitation of electrons is a small effect, a fact 
that holds if T <C Tp, where Tp is the Fermi temperature. 
For the densities of interest, those at which molecules dis- 
sociate in the fluid, Tp > 200, OOOK. The quahty of the 
BOA applied to high-pressure hydrogen fluid has been 
tested, and found to be good, by comparing CEIMC to 
accurate RPIMC calculations, which contain no such ap- 
proximation, at temperatures of 5,000 and 10, OOOK. [l^ 
We therefore anticipate that the BOA will be even more 
accurate at lower temperatures. 

For classical nuclei in the canonical ensemble, one can 
generate the configuration-space probability distribution 
by proposing a move from configuration sflE^ to con- 
figuration S" with uniform probability and using the 
Metropolis acceptance probability 



A = min[l,exp(-/3A(S', S"))] , 



(1) 



where (3 is the inverse temperature, (kBT)~ , and 
A (S*, iS") the difference in BO energies of nuclear con- 
figurations S and S", computed with a T = OK QMC 
method. The bias resulting from the statistical noise of 
A is reduced with correlated sampling [lo| , and elimi- 
nated by using the penalty method 1^. 

For electronic sampling, we use either VMC or rep- 
tation quantum Monte Carlo (RQMC). RQMC [l3| is 
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a projector method that does not suffer from the mixed- 
estimator bias of DMC and allows efficient calculations of 
energy differences. To handle electron antisymmetry, we 
use the fixed-phase method [l^ which allows the modulus 
of the many-electron wavefunction to be fully optimized 
while the phase is held fixed. This gives an estimate of 
the ground state energy much lower than that of the vari- 
ational estimate (by typically 2500K/atom) but above 
that of the exact ground state energy. Analysis of re- 
lated electron systems suggest that the error of the fixed- 
phase energy will be between 10% and 30% of the VMC 
error, (2^ i.e., between 250K/atom and 750K/atom. We 
expect relative errors, between different proton configura- 
tions, to be even less.jlH We also employ twist-averaged 
boundary conditions (TABC)[l9| in the electronic cal- 
culation to greatly reduce finite-size effects, a problem 
that has caused substantial errors in earlier simulations 
of fluid hydrogen. [231 Well-converged energies and cor- 
relation functions are achieved using, depending on the 
density, between 108 and 500 twist angles. 

In the CEIMC approach, one is free to choose any trial 
wavefunction for the electronic ground-state. An impor- 
tant consideration is the computer time required per step, 
since the wavefunction must be calculated for many nu- 
clear configurations. When studying dense liquid hydro- 
gen, we require a general wavefunction that is equally 
accurate for both molecular and non-molecular configu- 
rations. Therefore, we employ a fast band-structure cal- 
culation with an effective one-electron potential designed 
to produce accurate single-particle orbitals for a Slater- 
Jastrow type wavefunction. Within the fixed-phase ap- 
proach, it is only the orbitals that affect the systematic 
bias. One such set of orbitals would be those computed 
within Kohn-Sham theory. However, full self-consistency 
for each proton displacement would be too expensive 
for generating a large number of nuclear configuration- 
space samples. Consequently, we use a bare electron- 
proton potential for the effective single-particle poten- 
tial. Further screening and correlation effects are intro- 
duced with a Slater- Jastrow wavefunction; the Jastrow 
factor is from the RPA pseudopotential, including the 
one-body (electron- proton) term.[l^ Such an approach 
is surprisingly good, as demonstrated by Figure [T] which 
shows VMC and RQMC total energies for five differ- 
ent frozen nuclear configurations. This trial function 
is comparable in quality to one using Kohn-Sham LDA 
orbitals in the Slater determinant for all configurations 
and densities tested, while being faster to generate. It 
is also more transferable than both localized molecular 
orbitals Q (which cannot be used for non- molecular sys- 
tems) and analytic backflowjSsj (which is highly accurate 
only for high-density metallic systems). 

We investigate the atomic-molecular transition in liq- 
uid hydrogen at fixed density (1.35 < rs < 1.55 where 
47rrf/3 = l/rie) and temperature (T = 2000K and 
1500K). The pressure, estimated using the virial the- 




Configuration Configuration 



FIG. 1; VMC and RQMC total energies for five different crys- 
tal configurations at two different densities (left graphs V — 
9.42a.u./atom or Ts = 1.31; right graphs V — 33.51a.u./atom 
or Ts = 2.0) using a number of different Slater-Jastrow 
wavefunctions. Configuration 1 is molecular, 2-5 are non- 
molecular. Owing to the variational principle, a lower energy 
implies a better wavefunction. DFT-LDA and bare electron- 
proton bands (see text) provide the best and most transferable 
orbitals for the Slater determinant. Gaussian for configura- 
tion 1 refers to localized molecular orbitals. 

orem, an approach that is accurate with RQMC due to 
the lack of a mixed-estimator bias, lies between 135GPa 
and 290GPa for this range of densities. Using a classical 
Monte Carlo simulation with an empirical potential. (loj 
the system is prepared either in a purely molecular or a 
purely atomic fluid initial configuration. During the sub- 
sequent simulation the system evolves to its equilibrium 
state, subject to overcoming any free-energy barriers dur- 
ing the simulation. A typical simulation has 14,000 ionic 
moves with a step size of 0. 006-0. 016Bohr. We collect 
statistics along the sequence of ionic states, particularly 
the proton-proton correlation functions which give in- 
sight into the state of the liquid through a distinctive 
peak at rpp 1.4Bohr when molecules are present. Fig- 
ure m shows a set of proton-proton correlation functions 
for simulations at several densities prepared with either 
a molecular or atomic fluid as the initial state. 

For a quantitative analysis, a method for estimat- 
ing the average number of molecules from the pair- 
correlation function at each phase point is required. 
A rough estimate would involve integrating the pair- 
correlation function up to the first minimum, but such 
an approach does not take into account the "baseline" 
caused by nearby molecules and unbound atoms; see 
Fig. O To remove the baseline, we fit each of the pair- 
correlation functions to the functional form: 

g (r) = Xg^ (r; {a}) + (1 - X) ffat (r; {7}) , (2) 

where {a, 7} are fitting parameters and g^ is a Gaussian 
centered on 1.3-1.5Bohr. gat is a fit to the next-nearest 
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FIG. 2: Proton-proton correlation functions for CEIMC sim- 
ulations in the canonical ensemble with 32 atoms at T = 
2000K. The BO energies are computed with VMC (left pane) 
or RQMC (right pane) and the wavefunction is a Slater- 
Jastrow type with bare electron-proton bands (see text). Grey 
lines are simulations initially prepared with a molecular fluid 
and black with an atomic fluid. Hysteresis is evident. Pat and 
Pm are the computed pressures for simulations prepared with 
an atomic and a molecular fluid respectively. All correlation 
functions are plotted to the same scale. 



neighbor distribution. We find that the quality of the fit 
is insensitive to the choice of gat ; the molecular fraction 
varies by no more than 3% for different choices. The 
molecular order parameter, A, the fraction of protons that 
are bound into a molecule, is plotted against pressure in 
Fig. Elfor T = 2000K. 

We find that the CEIMC simulations using VMC en- 
ergy differences (left pane Fig. [21 and Fig. [3]) yield an ir- 
reversible transition with hysteresis. This implies a free- 
energy barrier that is difficult to overcome during the 
simulations, hence a metastable state is obtained. The 
large width of the hysteresis loop is indicative of a transi- 
tion that is weakly first-order. We therefore expect that 
T = 2000K is close to the VMC critical point. 

In contrast, the much more accurate RQMC simula- 
tions (right pane Fig. ^ and Fig. [3]) yield a reversible and 
continuous molecular dissociation upon increasing pres- 
sure. Both molecular and atomic states are mechanically 
unstable in the range of densities simulated, and the sys- 
tem quickly evolves, with no discernible free-energy bar- 
rier, to a stable part-molecular, part-atomic state. This 
behavior implies the lack of an underlying first-order 
transition, a behavior which would not change as the 
thermodynamic limit is approached. The likely reasons 
for the change in character of the crossover is that the BO 
energy surface of the assumed Slater-Jastrow trial func- 
tion is inaccurate under conditions at which molecules are 
short-lived transient entities. As discussed above, these 
deficiencies are corrected by the imaginary-time projec- 
tion of RQMC, if they are due to errors in the modulus 
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FIG. 3: Molecular order parameter (A in Eq. [2} evalu- 
ated for simulation results from Fig. [S] Two curves are 
obtained for each method from different initial conditions 
(see text). VMC (dashed lines) simulations indicate an ir- 
reversible, weakly first-order phase transition for molecu- 
lar dissociation in the fluid with increasing/decreasing pres- 
sure. This picture is compatible with conclusions from CPMD 
DFT-LDA. Accurate RQMC (solid lines) simulations flnd a 
reversible crossover. 

of the trial wavefunction and not its phase. 

We have assessed the finite-size errors in the cell vol- 
ume by repeating a selection of simulations with 54 
atoms. The proton-proton correlation functions in Figure 
|4] show the comparison between 32- and 54-atom simula- 
tions for VMC and RQMC. For VMC finite-size errors are 
clearly large, consistent with the presence of a free-energy 
barrier, while for RQMC the errors are small, indicative 
of a continuous crossover between molecular and atomic 
states. 

We have tested the nature of molecular dissociation 
at a lower temperature: T = 1500K. Figure \5\ shows 
pair-correlation functions. Hysteresis is again evident 
for VMC simulations, although over a shorter pressure- 
range, indicating a larger free-energy barrier. Again, 
RQMC simulations do not show a first-order transition, 
mirroring the behavior at T — 2000K. 

In conclusion, this study demonstrates that the 
CEIMC approach with a fast and accurate trial wave- 
function can be used to determine the nature of the 
pressure-driven molecular to atomic crossover in hydro- 
gen for temperatures on the order of lOOOK. We find 
that the improvements in interatomic interactions ob- 
tained when using RQMC to simulate the electronic 
subsystem lead to no first-order phase transition, a re- 
sult of fundamentally different nature to that predicted 
when using VMC. Our most accurate simulations, those 
employing RQMC, provide no evidence for a first-order 
atomic-molecular transition in the liquid phase at either 
T = 1500K or T = 2000K, showing instead a continu- 
ous crossover; this is in contrast to the first-order transi- 
tion found with CPMD. The remaining uncontrolled ap- 
proximations in our simulations are the finite-size error 
(demonstrated to be small for RQMC), the fixed-phase 
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FIG. 4: Demonstration of finite-size errors in VMC pair- 
correlation functions (upper pane). Sofid fines are 32-atoni 
simulations, dashed are 54 atoms. Grey lines are simulations 
started from a molecular fluid and black are from an atomic 
fluid. Quoted pressures are for simulations started from a 
molecular fluid. RQMC (lower pane) simulations show small 
flnite-size errors. 
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FIG. 5: Proton-proton correlation functions for simulations 
at r = 1500K. Black lines are simulations prepared with 
an atomic fluid, grey are prepared with a molecular fluid. 
The simulations demonstrate an irreversible transition with 
hysteresis for VMC. Left pane: VMC; Right pane: RQMC. 

approximation in RQMC, the adiabatic approximation 
for decoupling nuclear and electronic time-scales (tested 
m ref. and the use of classical statistics for simu- 

lating nuclei. Path integral calculations are in progress 
for studying the effect of nuclear quantum statistics on 
dense liquid hydrogen. We anticipate that nuclear quan- 



tum effects will destabilize molecules at a lower pressure 
leaving the crossover character unchanged. 
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